#!/bin/ksh

###################      AMB EL FONDO BLANC     ########################
####  		 REPRESENTACIO DE TRES FIGURES  	 ###############
########################################################################

DIR=`pwd`
##### input ########
#ftopo=/usr/users/ivone/Alps/topo/topo_filtrat.grd
ftopo=/home/ivone/Alps/topo/topo_EUR.grd
ffile1=e_sedsL_Tm_vis_Q_epeff_epzz.xy
fcorba1=corba.tmp
##### output ########
FILEPS=esLQ_ll.ps
rm $FILEPS

graficsth

Tcolor=1

awk '{print $0}' limits.d | read x1 x2 y1 y2 xtall ytall x2vel \
	y2vel VLLEG strlleg zx zy

awk '{print 0,0}' limits.d > file_xy.tmp
awk '{print $2*1e3+$1*1e3,$4*1e3+$3*1e3}' limits.d >> file_xy.tmp
awk '{if($1==">" || $1=="#") {print $0} else {print $1*1e3,$2*1e3}}' $fcorba1 >> file_xy.tmp
wc $fcorba1 | read npunts1 a b c
echo $npunts1+2 | bc -l | read nrow1
#wc $fcorba2 | read npunts2 a b c
#echo $nrow1+$npunts2 | bc -l | read nrow2
xytoll_reals.job	## input: file_xy.tmp,   output: file_ll.tmp
awk '{if(NR==1){print $0} }' file_ll.tmp | read lon1 lat1
awk '{if(NR==2){print $0} }' file_ll.tmp | read lon2 lat2
awk '{if(NR>2){print $0} }' file_ll.tmp > filecorba1.tmp

regiolon=$lon1/$lon2/$lat1/$lat2
regiol=2/23.4/44/53.4
#regiol=2/18.8/44/53.4
regiox=$x1/$x2/$y1/$y2

echo "--------------------------------------------------------------"
echo "    Regio real x, y:             " $regiox
echo "    Regio real longitut, latitut:" $regiolon
echo "    Regio figura lon, lat:       " $regiol
echo "--------------------------------------------------------------"
echo ' es coherent?'
read aaa
Dx_l=0.2
############################################################################
##############  PRIMERA - elevacio ##########################################

awk '{if($1!="#"){print $1*1e3,$2*1e3,$3} }' $ffile1 > file_xy.tmp
xytoll_reals.job	## input: file_xy.tmp,   output: file_ll.tmp

blockmedian file_ll.tmp -R$regiol -I$Dx_l -V > file_bloc.tmp
surface file_bloc.tmp -R -Gfile_grd.tmp -I$Dx_l -V
if [ $Tcolor -eq 1 ]
then
cat <<END>file_palette_cpt.tmp
-3     0  100  220    -2     0  150  240	
-2     0  150  220    -1     0  200  240	
-1     0  200  240     -0.500     0  220  255	
-0.5     0  220  255        0   120  255  255	
0    75  180  155      0.100    75  200  125	
0.100    75  200  125      0.200    75  235   75	
0.200    75  235   75      0.500   125  255   75	
0.500   150  255   75     1   175  255   75	
1   200  255   75     1.5   200  200   75  
1.5   180  160  100     2   180  160  100	
2   255  255  255     3.5   255  255  255	
END
else
cat <<END>file_palette_cpt.tmp
0  	10      10      10      0.4	10      10      10
0.4 	60      60      60      0.8	60      60      60
0.8	100     100     100     1.2	100     100     100
1.2	150     150     150     1.6	150     150     150
1.6	200     200     200     2.0	200     200     200
2.0	230     230     230     2.1	230     230     230
2.1	255     255     255     3.5	255     255     255
END
fi 
grdmath $ftopo 1000 DIV = file_topo_grd_km.tmp
pstext titol.tmp -P -R$regiox -X1 -Y25.1 -Jx.005 -N -K -V > $FILEPS
psbasemap -Y-4.3 -X1 -R$regiol -Jm0.42 -Ba4f2/a2f1WNes -K -O -V >> $FILEPS
##grdsample file_grd.tmp -Gfile_grd_sample.tmp -N101/101 -V 
file_grd=file_grd.tmp
#grd2cpt $file_grd > file_palette_cpt.tmp
grdimage file_topo_grd_km.tmp -Jm -R -O -K -Cfile_palette_cpt.tmp -V >> $FILEPS
pscoast -R -Jm -Dl -A4 -W4/255 -K -O -V >> $FILEPS
psxy filecorba1.tmp -R -O -K -M -Jm -W5/0/0/255t35_5_0_0:5 -V >> $FILEPS
grdcontour $file_grd -Ba -C0.3 -A0.6f6 -G2.5/8 -W3/0 -O -K -Jm -R -V >> $FILEPS
psscale -Cfile_palette_cpt.tmp -L -D13/2.7/5/.3 -B:."elevation (km)": -O -K -V >> $FILEPS
pstext -R -Jm -O -K -N -V <<END>> $FILEPS
25 54.5 9 0 4 0  $DIR
34 52.5 12 0 4 2  elevation (km)
END
#grdcontour $ftopo -C0.5 -A1f6 -G2.5/8 -W1/0 -O -K -Jm -R -V >> $FILEPS

rm file_*.tmp
#########################################################################
#####################  SEGONA -  gruix cortical ###############################
awk '{if($1!="#"){print ($1*1e3,$2*1e3,$5)} }' $ffile1 > file_xy.tmp
xytoll_reals.job	## input: file_xy.tmp,   output: file_ll.tmp

blockmedian file_ll.tmp -R$regiol -I$Dx_l -V > file_bloc.tmp
surface file_bloc.tmp -R -Gfile_grd.tmp -I$Dx_l -V
if [ $Tcolor -eq 1 ]
then
cat <<END>file_palette_cpt.tmp
20	255	0	255	24	255	0	255
24	113	0	255	26	113	0	255
26	0	28	255	28	0	28	255
28	0	170	255	30	0	170	255
30	0	255	199	32	0	255	199
32	0	255	56	34	0	255	56
34	85	255	0	36	85	255	0
36	227	255	0	38	227	255	0
38	255	142	0	40	255	142	0
40	255	0	0	60	255	0	0
END
else
cat <<END>file_palette_cpt.tmp
20  	0       0       0       24	0       0       0
24 	28      28      28      26	28      28      28
26 	57      57      57      28	57      57      57
28 	85      85      85      30	85      85      85
30	113     113     113     32	113     113     113
32	142     142     142     34	142     142     142
34	170     170     170     36	170     170     170
36	200     200     200     38	200     200     200
38	227     227     227     40	227     227     227
40	255     255     255     60	255     255     255
END
fi
psbasemap -Y-6.4 -R -Jm -Ba4f2/a2f1Wnes -K -O -V >> $FILEPS
file_grd=file_grd.tmp
grdimage $file_grd -Jm -R -O -K -Cfile_palette_cpt.tmp -V >> $FILEPS
pscoast -R -Jm -Dl -A4 -W4/255 -K -O -V >> $FILEPS
psxy filecorba1.tmp -R -O -K -M -Jm -W5/0/0/255t35_5_0_0:5 -V >> $FILEPS
grdcontour $file_grd -Ba -C2.5 -A5f6 -G2.5/8 -O -K -Jm -R -V >> $FILEPS
psscale -Cfile_palette_cpt.tmp -L -D13/2.7/5/.3 -B:."crustal thickness (km)  ": -O -K -V >> $FILEPS
pstext -R -Jm -O -K -N -V <<END>> $FILEPS
34 52.5 12 0 4 2 crustal thickness (km) 
END

rm file_*.tmp 

#########################################################################
################  TERCERA  - gruix litosferic ##########################
if [ $Tcolor -eq 1 ]
then
cat <<END>file_palette_cpt.tmp
70	255	0	255	80	255	0	255
80	113	0	255	85	113	0	255
85	0	28	255	90	0	28	255
90	0	170	255	95	0	170	255
95	0	255	199	100	0	255	199
100	0	255	56	105	0	255	56
105	85	255	0	110	85	255	0
110	227	255	0	115	227	255	0
115	255	142	0	120	255	142	0
120	255	0	0	150	255	0	0
END
else
cat <<END>file_palette_cpt.tmp
70  	0       0       0       75	0       0       0
75 	28      28      28      80	28      28      28
80 	57      57      57      85	57      57      57
85 	85      85      85      90	85      85      85
90	113     113     113     95	113     113     113
95	142     142     142     100	142     142     142
100	170     170     170     105	170     170     170
105	200     200     200     110	200     200     200
110	227     227     227     115	227     227     227
115	255     255     255     150	255     255     255
END
fi

awk '{if($1!="#"){print ($1*1e3,$2*1e3,$6)} }' $ffile1 > file_xy.tmp 
xytoll_reals.job	## input: file_xy.tmp,   output: file_ll.tmp

blockmedian file_ll.tmp -R$regiol -I$Dx_l -V > file_bloc.tmp
surface file_bloc.tmp -R -Gfile_grd.tmp -I$Dx_l -V

psbasemap -R -Y-6.4 -Jm -Ba4f2/a2f1Wnes -K -O -V >> $FILEPS
#grdfilter file_grd.tmp -D0 -Fb200 -Gfile_grd_filt.tmp -V
#grdsample litosfera_grd.tmp -Gfile_grd_sample.tmp -N100/100 -V -R
file_grd=file_grd.tmp
grdimage $file_grd -Jm -R -O -K -Cfile_palette_cpt.tmp -V >> $FILEPS
pscoast -R -Jm -Dl -A4 -W4/255 -K -O -V >> $FILEPS
psxy filecorba1.tmp -R -O -K -M -Jm -W5/0/0/255t35_5_0_0:5 -V >> $FILEPS
grdcontour $file_grd -Ba -C10 -A20f6 -G2.5/8 -O -K -Jm -R -V >> $FILEPS
psscale -Cfile_palette_cpt.tmp -L -D13/2.7/5/.3 -B:."lithospheric thickness (km)": -O -K -V >> $FILEPS
pstext -R -Jm -O -K -N -V <<END>> $FILEPS
34 52.5 12 0 4 2 lithospheric thickness (km)
END

grdsample file_grd.tmp -Gfile_2grd.tmp -I44.108/33.358 -R -V
grd2xyz file_2grd.tmp > file_L.tmp
sort -k 2,2n -k 1,1n -o fileL_ord.tmp file_L.tmp
rm file_*.tmp

#########################################################################
################  QUARTA - Flux de calor superficial ##########################
awk '{if($1!="#"){print ($1*1e3,$2*1e3,$9)} }' $ffile1 > file_xy.tmp
xytoll_reals.job	## input: file_xy.tmp,   output: file_ll.tmp

blockmedian file_ll.tmp -R$regiol -I$Dx_l -V > file_bloc.tmp
surface file_bloc.tmp -R -Gfile_grd.tmp -I$Dx_l -V
if [ $Tcolor -eq 1 ]
then
cat <<END>file_palette_cpt.tmp
50	255	0	255	55	255	0	255
55	113	0	255	60	113	0	255
60	0	28	255	65	0	28	255
65	0	170	255	70	0	170	255
70	0	255	199	75	0	255	199
75	0	255	56	80	0	255	56
80	85	255	0	85	85	255	0
85	227	255	0	90	227	255	0
90	255	142	0	95	255	142	0
95	255	0	0	100	255	0	0
END
else
cat <<END>file_palette_cpt.tmp
50  	0       0       0       55	0       0       0
55 	28      28      28      60	28      28      28
60 	57      57      57      65	57      57      57
65 	85      85      85      70	85      85      85
70	113     113     113     75	113     113     113
75	142     142     142     80	142     142     142
80	170     170     170     85	170     170     170
85	200     200     200     90	200     200     200
90	227     227     227     95	227     227     227
95	255     255     255     100	255     255     255
END
fi

psbasemap -R -Y-6.4 -Jm -Ba4f2/a2f1WSen -K -O -V >> $FILEPS
#grdfilter file_grd.tmp -D0 -Fb200 -Gfile_grd_filt.tmp -V
file_grd=file_grd.tmp
grdimage $file_grd -Jm -R -O -K -Cfile_palette_cpt.tmp -V >> $FILEPS
pscoast -R -Jm -Dl -A4 -W4/255 -K -O -V >> $FILEPS
psxy filecorba1.tmp -R -O -K -M -Jm -W5/0/0/255t35_5_0_0:5 -V >> $FILEPS
grdcontour $file_grd -Ba -C5 -A10f6 -G2.5/8 -O -K -Jm -R -V >> $FILEPS
psscale -Cfile_palette_cpt.tmp -L -D13/2.7/5/.3 -B:." ": -O -K -V >> $FILEPS
pstext -R -Jm -O -N -V <<END>> $FILEPS
34 52.5 12 0 4 2 Surface heat flow (mW/m@+2@+)
END

rm file*.tmp 
rm corba.tmp titol.tmp TITOL.tmp punts_falla_xy.res velocity.xy limits.d
rm e_sedsL_Tm_vis_Q_epeff_epzz.xy epunt12zz.xy escalar_strainrate.xy principal_stress.xy

gv -portrait -a4 -magstep -3 $FILEPS &
